/* Definitions for arithmetic types and functions.

   This file is part of khipu.

   khipu is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published by
   the Free Software Foundation; either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.  */

#include <cmath>
#include <climits>
#include "arith.hpp"
#include "floatp.hpp"

KP_DECLS_BEGIN

// Assumes 'interp' is bound.
#define kp_allocal(n)   \
  (limb_t *)KP_TALLOC (interp, (n) * sizeof (limb_t))

int uicmpn (const limb_t *xp, const limb_t *yp, int n)
{
  for (; n > 0; --n)
    if (xp[n - 1] > yp[n - 1])
      return (+1);
    else if (xp[n - 1] < yp[n - 1])
      return (-1);

  return (0);
}

int uicmp (const limb_t *xp, int xl,
  const limb_t *yp, int yl)
{
  int t = xl - yl;
  return (t ? t : uicmpn (xp, yp, xl));
}

limb_t uiadd1 (limb_t *dst,
  const limb_t *src, int n, limb_t val)
{
  for (int i = 0; i < n; ++i)
    {
      limb_t prev = src[i];
      dst[i] = prev + val;
      if (!(val = dst[i] < prev))
        {
          if (dst != src)
            memcpy (&dst[i + 1], &src[i + 1],
              (n - i - 1) * sizeof (limb_t));
          break;
        }
    }

  return (val);
}

int uiaddn (limb_t *dst, const limb_t *xp,
  const limb_t *yp, int n)
{
  limb_t cy = 0;

  for (int i = 0; i < n; ++i)
    {
      dst[i] = xp[i] + yp[i] + cy;
      cy = !cy ? dst[i] < yp[i] : dst[i] <= yp[i];
    }

  return ((int)cy);
}

int uiadd (limb_t *dst, const limb_t *xp,
  int xl, const limb_t *yp, int yl)
{
  if (xl < yl)
    {
      swap (xp, yp);
      swap (xl, yl);
    }

  if (uiaddn (dst, xp, yp, yl) != 0)
    {
      do ;
      while (++dst[yl++] == 0);

      return (yl > xl);
    }

  if (dst != xp)
    memcpy (&dst[yl], &xp[yl],
      (xl - yl) * sizeof (limb_t));

  return (0);
}

limb_t uisub1 (limb_t *dst,
  const limb_t *src, int n, limb_t val)
{
  for (int i = 0; i < n; ++i)
    {
      int tmp = src[i] < val;
      dst[i] = src[i] - val;
      if (!(val = tmp))
        {
          if (dst != src)
            memcpy (&dst[i + 1], &src[i + 1],
              (n - i - 1) * sizeof (limb_t));

          val = i == n - 1;
          break;
        }
    }

  return (val);
}

int uisubn (limb_t *dst, const limb_t *xp,
  const limb_t *yp, int n)
{
  limb_t cy = 0;

  for (int i = 0; i < n; ++i)
    {
      limb_t prev = xp[i];
      dst[i] = prev - yp[i] - cy;
      cy = !cy ? prev < yp[i] : prev <= yp[i];
    }

  return ((int)cy);
}

int uisub (limb_t *dst, const limb_t *xp,
  int xl, const limb_t *yp, int yl)
{
  if (uisubn (dst, xp, yp, yl))
    {
      limb_t v;

      do
        {
          if (yl >= xl)
            return (1);

          v = xp[yl];
          dst[yl++] = v - 1;
        }
      while (v == 0);
    }

  if (dst != xp)
    memcpy (dst + yl, xp + yl, (xl - yl) * sizeof (limb_t));

  return (0);
}

int uilsh (limb_t *dst, const limb_t *src, int len, int shift)
{
  uival acc = { 0 };
  int i;

  for (i = 0; i < len; ++i)
    {
      acc.qv += (twolimb_t)src[i] << shift;
      dst[i] = acc.limbs.lo;
      acc.cuthi ();
    }

  return ((dst[i] = acc.limbs.lo) != 0);
}

int uimul2exp (limb_t *dst, const limb_t *src, int len, int shift)
{
  int cnt = shift / LIMB_BITS;
  int ret = len + cnt;

  if ((shift %= LIMB_BITS) != 0)
    ret += uilsh (dst + cnt, src, len, shift);
  else
    memmove (dst + cnt, src, len * sizeof (*dst));

  memset (dst, 0, cnt * sizeof (*dst));
  return (ret);
}

int uirsh (limb_t *dst, const limb_t *src, int len, int shift)
{
  limb_t *outp = dst + len;
  const limb_t *inp = src + len;
  uival acc = { 0 };

  while (inp != src)
    {
#ifdef KP_ARCH_WIDE
      acc.qv <<= LIMB_BITS;
#else
      acc.limbs.hi = acc.limbs.lo, acc.limbs.lo = 0;
#endif
      acc.qv += (twolimb_t)*--inp << (LIMB_BITS - shift);
      *--outp = acc.limbs.hi;
    }

  return (dst[len - 1] != 0);
}

int uidiv2exp (limb_t *dst, const limb_t *src, int len, int shift)
{
  int cnt = shift / LIMB_BITS;
  int ret = len - cnt;

  if (ret <= 0)
    ret = 0;
  else if ((shift %= LIMB_BITS) == 0)
    memmove (dst, src + cnt, ret * sizeof (*dst));
  else if (!uirsh (dst, src + cnt, ret, shift))
    --ret;

  return (ret);
}

int uimul1 (limb_t *dst, const limb_t *src, int len, limb_t val)
{
  limb_t prev = 0;
  int i;

  for (i = 0; i < len; ++i)
    {
      uival cy;

      cy.qmuladd (val, src[i], prev);
      dst[i] = cy.limbs.lo;
      prev = cy.limbs.hi;
    }

  return ((dst[i] = prev) != 0);
}

static const int KARATSUBA_THRESHOLD = 25;

limb_t uimul (interpreter *interp, limb_t *ret,
  const limb_t *xp, int xl, const limb_t *yp, int yl)
{
  if (kp_likely (yl <= KARATSUBA_THRESHOLD || xl <= KARATSUBA_THRESHOLD))
    { // Basic multiplication.
      for (int i = 0; i < xl; ++i)
        {
          uival cy = { 0 };
          for (int j = 0; j < yl; ++j)
            {
              uival tmp;

              tmp.qmuladd (xp[i], yp[j], ret[i + j]);
              cy.qv += tmp.qv;
              ret[i + j] = cy.limbs.lo;
              cy.cuthi ();
            }

          ret[i + yl] = cy.limbs.lo;
        }

      return (ret[xl + yl - 1]);
    }

  // Very large numbers - Use Karatsuba algorithm.
  int shift = (xl + yl) * (LIMB_BITS >> 2);
  int word = shift / LIMB_BITS;
  int bytes = (1 << (shift % LIMB_BITS)) - 1;

  KP_TMARK (interp);
  limb_t *x0p = kp_allocal (word * 2 + xl + yl);

  // x0 = x & ((1 << shift) - 1) 
  memcpy (x0p, xp, word * sizeof (limb_t));
  x0p[word - 1] &= bytes;

  // x1 = x >> shift */
  limb_t *x1p = x0p + word;
  int x1l = uidiv2exp (x1p, xp, xl, shift);

  // y0 = y & ((1 << shift) - 1)
  limb_t *y0p = x1p + xl;
  memcpy (y0p, yp, word * sizeof (limb_t));
  y0p[word - 1] &= bytes;

  // y1 = y >> shift
  limb_t *y1p = y0p + word;
  int y1l = uidiv2exp (y1p, yp, yl, shift);

  limb_t *z0p = kp_allocal (word * 2 + x1l + y1l);
  memset (z0p, 0, (word * 2 + x1l + y1l) * sizeof (limb_t));
  limb_t *z2p = z0p + word * 2;

  // z0 = x0 * y0, z2 = x1 * y1
  int z2l = x1l + y1l, z0l = word * 2;
  z0l -= uimul (interp, z0p, x0p, word, y0p, word) == 0;
  z2l -= uimul (interp, z2p, x1p, x1l, y1p, y1l) == 0;

  // y1 += y0, x1 += x0 
  int rl = uiadd (y1p, y1p, y1l, y0p, word);
  y1l += rl;
  rl = uiadd (x1p, x1p, x1l, x0p, word);
  x1l += rl;

  // z1 = (x0 + x1) * (y0 + y1)
  int z1l = x1l + y1l;
  limb_t *z1p = kp_allocal (z1l);
  memset (z1p, 0, z1l * sizeof (limb_t));
  z1l -= uimul (interp, z1p, x1p, x1l, y1p, y1l) == 0;

  // z1 -= (z2 + z0)
  z1l -= uisub (z1p, z1p, z1l, z2p, z2l);
  z1l -= uisub (z1p, z1p, z1l, z0p, z0l);

  // ret = (z1 << (shift << 1)) + (z1 << shift) + z0
  rl  = uimul2exp (ret, z2p, z2l, shift);
  rl += uiadd (ret, ret, rl, z1p, z1l);
  rl  = uimul2exp (ret, ret, rl, shift);
  uiadd (ret, ret, rl, z0p, z0l);

  return (ret[xl + yl - 1]);
}

limb_t uidivrem1 (limb_t *dst,
  const limb_t *src, int len, limb_t val)
{
  limb_t *dp = dst + len;
  const limb_t *sp = src + len;
  limb_t rem = 0;

  for (int i = 0; i < len; ++i)
    {
      uival tmp;
#ifdef KP_ARCH_WIDE
      tmp.qv = ((twolimb_t)rem << LIMB_BITS) | *--sp;
#else
      tmp.limbs.hi = rem, tmp.limbs.lo = *--sp;
#endif
      *--dp = (limb_t)(tmp.qv / val);
      rem = (limb_t)(tmp.qv % val);
    }

  return (rem);
}

int uidiv (interpreter *interp,
  const limb_t *xp, int xl,
  const limb_t *yp, int yl,
  limb_t *roomp, bool rem_p)
{
  limb_t *r = roomp, *q = r + yl;
  int rl = yl, ql = xl - yl + 1;
  limb_t s = uiclz (yp[yl - 1]);
  KP_TMARK (interp);

  if (s > 0)
    { // Left shift Y so that the MSB is set. Do the same for X.
      limb_t *dp = kp_allocal (yl);
      uilsh (dp, yp, yl, s);
      uilsh (r, xp, xl, s);
      yp = dp;
    }
  else
    memcpy (r, xp, xl * sizeof (limb_t));

  limb_t ymsd = yp[yl - 1];
  limb_t y2msd = yp[yl - 2];
  twolimb_t ymsdd = ((twolimb_t)ymsd << LIMB_BITS) | y2msd;

  // Division loop.
  for (int i = ql - 1; i >= 0; --i)
    {
      limb_t qstar, c1;

      if (q[i] < ymsd)
        { // No overflow - Compute initial q*.
          twolimb_t num = ((twolimb_t)q[i] << LIMB_BITS) | q[i - 1];
          qstar = (limb_t)(num / ymsd);
          c1 = (limb_t)(num % ymsd);
        }
      else
        { // Overflow - Adjust q*.
          qstar = ~(limb_t)0;
          if (q[i] > ymsd || (c1 = ymsd + q[i - 1]) < ymsd)
            goto sub;
        }

      {
        twolimb_t c2 = ((twolimb_t)c1 << LIMB_BITS) | q[i - 2];
        uival c3;

        c3.qmul (y2msd, qstar);
        // We know this can happen at most twice, so adjust q* accordingly.
        if (c3.qv > c2)
          qstar -= 1 + (c3.qv - c2 > ymsdd);
      }

      if (qstar > 0) sub:
        { // Subtract q* from the remainder.
          limb_t *dp = &r[i];
          uival cy = { 0 };

          for (int j = 0; j < yl; ++j)
            {
              uival tmp;

              tmp.qmuladd (qstar, yp[j], ~*dp);
              cy.qv += tmp.qv;
              *dp++ = ~tmp.limbs.lo;
              cy.cuthi ();
            }

          if (cy.limbs.lo > q[i])
            {
              limb_t cy2 = 0;
              dp -= yl, --qstar;
              for (int j = 0; j < yl; ++j)
                {
                  limb_t s1 = yp[j], s2 = *dp;
                  *dp++ = s1 + s2 + cy2;
                  cy2 = !cy2 ? s1 > ~s2 : s1 >= ~s2;
                }
            }
        }

      // Store q* in the quotient.
      q[i] = qstar;
    }

  if (rem_p && s > 0)
    // Right shift R back only if the remainder was requested.
    uirsh (r, r, rl, s);

  uitrim (q, ql);
  return (ql);
}

limb_t uisqr (interpreter *interp,
  limb_t *dst, const limb_t *src, int len)
{
  int i, imax = len;
  limb_t *tmp;
  KP_TMARK (interp);

  for (i = 1; i < imax; ++i)
    {
      uival k = { 0 };
      
      for (int j = 0; j < i; ++j)
        {
          uival tmp;

          tmp.qmuladd (src[j], src[j], dst[i + j]);
          k.qv += tmp.qv;
          dst[i + j] = k.limbs.lo;
          k.cuthi ();
        }
      
      dst[i + i] = k.limbs.lo;
    }

  len <<= 1;
  uitrim (dst, len);
  len += uilsh (dst, dst, len, 1);
  
  tmp = kp_allocal (imax + imax);
  for (i = 0; i < imax; ++i)
    {
      uival k;

      k.qmul (src[i], src[i]);
      tmp[2 * i + 0] = k.limbs.lo;
      tmp[2 * i + 1] = k.limbs.hi;
    }

  i <<= 1;
  uitrim (tmp, i);
  uiadd (dst, dst, len, tmp, i);
  return (dst[imax * 2 - 1]);
}

// For 2 <= Base <= 32.
static const int LOG2TAB[] =
{
  1, 0, 2, 0, 0, 0, 3, 0, 0, 0,
  0, 0, 0, 0, 4, 0, 0, 0, 0, 0,
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5 
};

static inline int uilog2 (int radix)
{
  return (LOG2TAB[radix - 2]);
}

// For 3 <= Base <= 36
/* `rbase' represents a number Beta, such that Beta = Base^N, N being
 * floor (log2b * LNBASE_BITS). In other words, Beta is the maximum power of
 * the radix that fits in a `limb_t'. It is used when generating digits:
 * Instead of dividing by the radix, we divide by Beta, obtaining exactly
 * N digits, and performing less divisions. */
static const struct
{
  limb_t rbase;
  double log2b;
} IBTAB[] =  {
  // The holes in the table belong to the powers of 2.
  { 0xcfd41b91, 0.6309297535714574 }, { 0 },
  { 0x48c27395, 0.4306765580733931 },
  { 0x81bf1000, 0.3868528072345416 },
  { 0x75db9c97, 0.3562071871080222 }, { 0 },
  { 0xcfd41b91, 0.3154648767857287 },
  { 0x3b9aca00, 0.3010299956639811 },
  { 0x8c8b6d2b, 0.2890648263178878 },
  { 0x19a10000, 0.2789429456511298 },
  { 0x309f1021, 0.2702381544273197 },
  { 0x57f6c100, 0.2626495350371936 },
  { 0x98c29b81, 0.2559580248098155 }, { 0 },
  { 0x18754571, 0.2446505421182260 },
  { 0x247dbc80, 0.2398124665681315 },
  { 0x3547667b, 0.2354089133666382 },
  { 0x4c4b4000, 0.2313782131597592 },
  { 0x6b5a6e1d, 0.2276702486969530 },
  { 0x94ace180, 0.2242438242175754 },
  { 0xcaf18367, 0.2210647294575037 },
  { 0x0b640000, 0.2181042919855316 },
  { 0x0e8d4a51, 0.2153382790366965 },
  { 0x1269ae40, 0.2127460535533632 },
  { 0x17179149, 0.2103099178571525 },
  { 0x1cb91000, 0.2080145976765095 },
  { 0x23744899, 0.2058468324604345 },
  { 0x2b73a840, 0.2037950470905062 },
  { 0x34e63b41, 0.2018490865820999 }, { 0 },
  { 0x4cfa3cc1, 0.1982398631705605 },
  { 0x5c13d840, 0.1965616322328226 },
  { 0x6d91b519, 0.1949590218937863 },
  { 0x81bf1000, 0.1934264036172708 }
};

int uibsize (int radix, limb_t msd, int len)
{
  if (pow2_p (radix))
    {
      int bits = len * LIMB_BITS - uiclz (msd);
      int bb = uilog2 (radix);
      return ((bits + bb - 1) / bb);
    }

  return (2 + len * (int)(LIMB_BITS * IBTAB[radix - 3].log2b));
}

int invbsize (int radix, int len)
{
  if (pow2_p (radix))
    {
      int bb = uilog2 (radix);
      return ((bb * (len - 1) + 1) / LIMB_BITS);
    }
  else
    return ((int)(len / (LIMB_BITS * IBTAB[radix - 3].log2b)));
}

static const char DIGUPR[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
static const char DIGLWR[] = "0123456789abcdefghijklmnopqrstuvwxyz";

int uitostr (interpreter *interp, char **dstp, int cap,
  const limb_t *src, int len, int radix)
{
  const char *digs = radix < 0 ? (radix = -radix, DIGUPR) : DIGLWR;
  char *dst = *dstp;

  if (pow2_p (radix))
    {
      limb_t n[] = { 0, src[len - 1] };
      int bpd = uilog2 (radix);
      int cnt = uiclz (n[1]);
      char *ptr = dst;
      int bits = LIMB_BITS * len - cnt;
      
      if (0 != (cnt = (int)(bits % bpd)))
        bits += bpd - cnt;
      
      int bpos = (int)(bits - (len - 1) * LIMB_BITS);
      for (int i = len - 1 ; ; )
        {
          for (bpos -= bpd; bpos >= 0; bpos -= bpd)
            *ptr++ = digs[(n[1] >> bpos) & ((1 << bpd) - 1)];
          
          if (--i < 0)
            break;
          
          n[0] = (n[1] << -bpos) & ((1 << bpd) - 1);
          n[1] = src[i];
          bpos += LIMB_BITS;
          *ptr++ = digs[n[0] | (n[1] >> bpos)];
        }
      
      return ((int)(ptr - dst));
    }
  else
    {
      char *dptr, *cptr = dst + cap - 1;
      limb_t ub = IBTAB[radix - 3].rbase;
      int base_exp = (int)((double)LIMB_BITS * IBTAB[radix - 3].log2b);

      // Input is clobbered when the radix isn't a power of 2.
      KP_TMARK (interp);
      limb_t *tmp = kp_allocal (len);
      memcpy (tmp, src, len * sizeof (limb_t));

      // Divide by radix^exp, in place.
      for (dptr = cptr; len > 0; )
        {
          limb_t rem = uidivrem1 (tmp, tmp, len, ub);
          
          for (int cnt = 0; cnt != base_exp; ++cnt, rem /= radix)
            *dptr-- = digs[rem % radix];
          
          if (!tmp[len - 1]) --len;
        }
      
      // Skip leading zeroes.
      while (*++dptr == '0') ;
      *dstp = dptr;
      return ((int)(cptr - dptr + 1));
    }
}

char* uitostr1 (char *dst, limb_t val, int radix)
{
  const char *digs = radix < 0 ? (radix = -radix, DIGUPR) : DIGLWR;

  do
    *--dst = digs[val % radix];
  while ((val /= radix) != 0);

  return (dst);
}

/* The digit value of a char C, equivalent to:
 * - for C in ['0' .. '9'] => C - '0'
 * - for C in ['a' .. 'z'] => C - 'a' + 10
 * - for C in ['A' .. 'Z'] => C - 'A' + 10 .*/
static const limb_t VINBASE[] =
{
  0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07,
  0x08, 0x09, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
  0xff, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10,
  0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18,
  0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0x20,
  0x21, 0x22, 0x23, 0xff, 0xff, 0xff, 0xff, 0xff,
  0xff, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0x10,
  0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18,
  0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0x20,
  0x21, 0x22, 0x23
};

static inline limb_t radixval (int c)
{
  return (VINBASE[(unsigned char)c -  '0']);
}

intptr_t strtoui1 (const char *src, int len, int radix)
{
  intptr_t ret = 0;

  for (int ix = 0; ix < len; ++ix)
    ret = ret * radix + radixval (src[ix]);

  return (ret);
}

int strtoui (limb_t *dst, const char *src, int len, int radix)
{
  int ret = 0;

  if (pow2_p (radix))
    {
      int next_bp = 0;
      limb_t digit = 0;
      int bpi = uilog2 (radix);
      
      for (const char *s = src + len - 1; s >= src; --s)
        {
          int inp = radixval (*s);

          digit |= ((limb_t)inp << next_bp);
          if ((next_bp += bpi) >= LIMB_BITS)
            {
              dst[ret++] = digit;
              next_bp -= LIMB_BITS;
              digit = (limb_t)inp >> (bpi - next_bp);
            }
        }
      
      if (digit != 0)
        dst[ret++] = digit;
    }
  else
    {
      limb_t limit = (~(limb_t)0) / radix;
      limb_t mul = 1, acc = 0;
      const char *end = src + len;
      
      for (++ret; src < end; ++src)
        {
          acc = acc * radix + radixval (*src);
		  
          if ((mul *= radix) >= limit)
            {
              limb_t cy = uimul1 (dst, dst, ret, mul);
              ret += cy;
              cy = uiadd1 (dst, dst, ret, acc);
              dst[ret] = cy, ret += cy;
              // Reset values.
              mul = 1, acc = 0;
            }
        }

      if (mul > 1)
        {
          limb_t cy = uimul1 (dst, dst, ret, mul);
          ret += cy;
          cy = uiadd1 (dst, dst, ret, acc);
          dst[ret] = cy, ret += cy;
        }
    }

  return (ret);
}

static int
uipow1_hi (interpreter *interp, limb_t *dst, limb_t *ignp,
  int base, int expo, int len, limb_t *tp)
{
  if (!expo)
    {
      *dst = 1, *ignp = 0;
      return (1);
    }

  int ign = 0, off = 0, rl = 1;
  int cnt = uiclz (expo);
  limb_t *rp = dst;

  *dst = base;
  for (int i = LIMB_BITS - cnt - 2; i >= 0; --i)
    {
      memset (tp, 0, rl * 2 * sizeof (limb_t));
      limb_t cy = uisqr (interp, tp, dst + off, rl);
      rl = 2 * rl - (cy == 0);
      ign <<= 1, off = 0;

      if (rl > len)
        {
          ign += (off = rl - len);
          rl = len;
        }
      
      swap (tp, dst);
      
      if (((expo >> i) & 1) != 0)
        {
          cy = uimul1 (dst, dst + off, rl, base);
          rl += cy, off = 0;
        }
    }

  if (rl > len)
    {
      ign += rl - len;
      dst += rl - len;
      rl = len;
    }

  if (rp != dst + off)
    memcpy (rp, dst + off, rl * sizeof (limb_t));

  *ignp = ign;
  return (rl);
}

static int
lfrounds (char *dst, int len, uint32_t base, int *expp)
{
  int off = 0;
  if (radixval (dst[len--]) * 2 >= base)
    for ( ; ; --len, ++off)
      {
        if (radixval (dst[len]++) != base - 1)
          break;
        else if (len == 0)
          {
            *dst = '1';
            ++*expp;
            break;
          }
      }

  return (off);
}

int lftostr (interpreter *interp, char *dst, int *expp,
  const limb_t *src, int len, int digs, int ibase)
{
  int expo = *expp;
  int needed_digs, nd, tl, exp_base;
  limb_t ign, base = abs (ibase);
  KP_TMARK (interp);

  int got, cap = digs + 2 * LIMB_BITS + 3;
  char *tmps = (char *)KP_TALLOC (interp, cap);

  if (pow2_p (base))
    {
      int bl2 = uilog2 (base);
      needed_digs = (digs * bl2) / LIMB_BITS + 2;
      nd = LIMB_BITS + abs (needed_digs - expo) / bl2;
    }
  else
    {
      double bl2 = IBTAB[base - 3].log2b;
      needed_digs = (int)(digs / bl2) / LIMB_BITS + 2;
      nd = (int)(LIMB_BITS * abs (needed_digs - expo) * bl2);
    }

  if (len > needed_digs)
    {
      src += len - needed_digs;
      len = needed_digs;
    }

  limb_t *pp = kp_allocal ((needed_digs + 1) * 4);
  limb_t *tp = &pp[(needed_digs + 1) * 2];
  int pl = uipow1_hi (interp, pp, &ign, base, nd, needed_digs, tp);

  if (expo <= needed_digs)
    { // Multiply by base^exp.
      int tl = len + pl, off = len - expo - ign;

      memset (tp, 0, tl * sizeof (limb_t));
      tl -= uimul (interp, tp, src, len, pp, pl) == 0;
      if (off < 0)
        {
          memcpy (tp - off, tp, tl * sizeof (limb_t));
          memset (tp, 0, -off * sizeof (limb_t));
          tl -= off, off = 0;
        }

      got = uitostr (interp, &tmps, cap, tp + off, tl - off, ibase);
      exp_base = got - nd;
    }
  else
    { // Divide by base^exp.
      int xl = expo - ign, off = xl - len;
      limb_t *xp;

      if (xl > len)
        { // Extend and zero-fill X.
          xp = kp_allocal (xl);
          memset (xp, 0, off * sizeof (limb_t));
          memcpy (xp + off, src, len * sizeof (limb_t));
        }
      else
        xp = (limb_t *)src, xl = len;

      tl = uidiv (interp, xp, xl, pp, pl, tp, false);
      got = uitostr (interp, &tmps, cap, tp + pl, tl, ibase);
      exp_base = got + nd;
    }

  // Round the result if more than requested digits were computed.
  if (got > digs)
    got = digs - lfrounds (tmps, digs, base, &exp_base);

  if (digs > got)
    digs = got;

  for (; digs > 0 && tmps[digs - 1] == '0'; --digs) ;
  memcpy (dst, tmps, digs);
  *expp = exp_base;
  return (digs);
}

static inline bool
frac_char_p (int ch, int radix)
{
  return (ch == '.' || ch == '@' ||
    (radix == 10 && (ch == 'e' || ch == 'E')));
}

int parse_num (interpreter *interp, const char *src,
  int len, num_info& info)
{
  int ix = 0;

  memset (&info, 0, sizeof (info));
  info.type = typecode::INT;
  info.radix = 10;

  if (src[ix] == '-')
    info.sign = 1, ++ix;
  else if (src[ix] == '+')
    ++ix;

  if (len == ix)
    return (-1);
  else if (src[ix] == '0')
    {
      if (len - ix <= 1)   // +/- 0
        return (0);
      else if (src[ix + 1] == 'x' || src[ix + 1] == 'X')
        info.radix = 16, ix += 2;
      else if (src[ix + 1] == 'b')
        info.radix = 2, ix += 2;
      else if (src[ix + 1] == 'r')
        {
          if (len - ix <= 4)
            return (-1);   // "+/-0r"

          info.radix = (src[ix + 2] - '0') * 10 + src[ix + 3] - '0';
          if (info.radix > 36)
            return (-1);

          ix += 3;
        }
      else if (src[ix + 1] != '.' && src[ix + 1] != 'e' && src[ix + 1] != 'E')
        info.radix = 8, ++ix;
    }

  for (; src[ix] == '0' && ix < len; ++ix) ;
  if (ix == len)
    return (0);   // Nothing but 0's.

  info.dec_start = ix;

#define CHR_OOB()   \
  src[ix] < '0' || (uint32_t)(src[ix] - '0') >= KP_NELEM (VINBASE) ||   \
      radixval (src[ix]) >= (uint32_t)info.radix

  for (; ix < len; ++ix)
    {
      if (frac_char_p (src[ix], info.radix))
        break;
      else if (CHR_OOB ())
        return (-1);
    }

  info.dec_end = ix;
  if (!frac_char_p (src[ix], info.radix))
    return (0);   // An integer value.

  info.type = typecode::FLOAT;
  info.got_dot = src[ix] == '.';
  for (ix += info.got_dot; ix < len; ++ix)
    {
      if (src[ix] == '@' || (info.radix == 10 &&
          (src[ix] == 'e' || src[ix] == 'E')))
        {
          info.frac_end = ix;
          break;
        }
      else if (CHR_OOB ())
        return (-1);
    }

  info.frac_end = ix;
  if (ix != len)
    { // We have an exponent.
      if (!ix || len - ++ix < 1)
        return (-1);
      else if (src[ix] == '-')
        info.expo_sign = 1, ++ix;
      else if (src[ix] == '+')
        ++ix;

      if (len - ix < 1)
        return (-1);

      for (; src[ix] == '0' && ix < len; ++ix) ;

      if (ix != len)
        {
          info.expo_start = ix;
          for (; ix < len; ++ix)
            if (CHR_OOB ())
              return (-1);
        }
    }

#undef CHR_OOB
  return (0);
}

static int
strtoui_lf_helper (interpreter *interp, const char *src,
  int slen, const num_info& info, limb_t *outp)
{
  KP_TMARK (interp);
  char *sp = (char *)KP_TALLOC (interp, slen), *up = sp;
  for (int i = info.dec_start; i < info.dec_end; ++i)
    *up++ = src[i];

  for (int i = info.dec_end + 1; i < info.frac_end; ++i)
    *up++ = src[i];

  return (strtoui (outp, sp, (int)(up - sp), info.radix));
}

int strtolf (interpreter *interp, const char *src, int slen,
  const num_info& info, int& out_expo, limb_t *outp)
{
  int expo = out_expo;
  int rn = 1 + invbsize (info.radix, abs (expo));
  int rl = strtoui_lf_helper (interp, src, slen, info, outp);

  if (expo == 0)
    return (out_expo = rl);
  else if (expo < 0)
    {
      int diff = rl - rn;
      int off = diff < 0 ? -diff + 1 : 1;
      int a_expo = -expo;

      memmove (outp + off, outp, rl * sizeof (*outp));
      memset (outp, 0, off * sizeof (*outp));
      rl += off;

      int base_exp = (int)(LIMB_BITS * IBTAB[info.radix - 3].log2b);
      limb_t ub = IBTAB[info.radix - 3].rbase;

      for (; a_expo > base_exp; a_expo -= base_exp)
        {
          uidivrem1 (outp, outp, rl, ub);
          if (outp[rl - 1] == 0)
            --rl;
        }

      expo = 1;
      for (int i = 0; i < a_expo; ++i, expo *= info.radix);
      uidivrem1 (outp, outp, rl, expo);
      if (outp[rl - 1] == 0)
        --diff, --rl;

      out_expo = diff + 1;
    }
  else
    {
      int a_expo = expo;
      int base_exp = (int)(LIMB_BITS * IBTAB[info.radix - 3].log2b);
      limb_t ub = IBTAB[info.radix - 3].rbase;

      for (; a_expo > base_exp; a_expo -= base_exp)
        if (uimul1 (outp, outp, rl, ub))
          ++rl;

      expo = 1;
      for (int i = 0; i < a_expo; ++i, expo *= info.radix);
      if (uimul1 (outp, outp, rl, expo))
        ++rl;

      out_expo = rl;
    }

  return (rl);
}

#ifndef __GNUC__

#ifdef _MSC_VER
  #include <intrin.h>
#  pragma intrinsic(_BitScanReverse)
#else

static int
uiones (limb_t x)
{
  x -= (x >> 1) & 0x55555555;
  x = ((x >> 2) & 0x33333333) + (x & 0x33333333);
  x = ((x >> 4) + x) & 0xF0F0F0F;
  x += x >> 8;
  x += x >> 16;
  return (x & 0x3F);
}

#endif   // _MSC_VER

unsigned int uiclz (limb_t val)
{
#ifdef _MSC_VER
  unsigned long ret;
  _BitScanReverse (&ret, (unsigned long)val);
#else
  val |= val >> 1;
  val |= val >> 2;
  val |= val >> 4;
  val |= val >> 8;
  val |= val >> 16;
  return (LIMB_BITS - uiones (val));
#endif
}

#endif   // !defined (__GNUC__)

#if FLT_RADIX == 2
#  define DBL_MANT_BIT   DBL_MANT_DIG
#elif FLT_RADIX == 4
#  define DBL_MANT_BIT   (DBL_MANT_DIG * 2)
#elif FLT_RADIX == 8
#  define DBL_MANT_BIT   (DBL_MANT_DIG * 3)
#elif FLT_RADIX == 16
#  define DBL_MANT_BIT   (DBL_MANT_DIG * 4)
#endif

int dbltoui (double src, limb_t *dst, int& len)
{
  int iexp;
  double y = frexp (src, &iexp);

  memset (dst, 0, LNDBL_SIZE * sizeof (*dst));

#if (DBL_MANT_BIT % KP_LIMB_BITS) != 0
#  if (DBL_MANT_BIT % KP_LIMB_BITS) > KP_LIMB_BITS / 2

  {
    y *= 1u << (DBL_MANT_BIT % (LIMB_BITS / 2));
    uint32_t hi = (uint32_t)y;
    y -= hi;

    y *= 1u << (LIMB_BITS / 2);
    uint32_t lo = (uint32_t)y;
    y -= lo;
    dst[DBL_MANT_BIT / LIMB_BITS] = (hi << (LIMB_BITS / 2)) | lo;
  }

#  else

  {
    y *= 1u << (DBL_MANT_BIT % LIMB_BITS);
    uint32_t d = (uint32_t)y;
    y -= d;
    dst[DBL_MANT_BIT / LIMB_BITS] = d;
  }

#  endif
#endif

  for (int i = DBL_MANT_BIT / LIMB_BITS; i > 0; )
    {
      y *= 4294967296.;
      uint32_t val = (uint32_t)y;
      dst[--i] = val;
      y -= val;
    }

  int rl = LNDBL_SIZE - 1;
  uitrim (dst, rl);
  int ret = rl;

  iexp -= DBL_MANT_BIT;   // The real exponent, in radix 2
  int a_exp = abs (iexp);

  /* Now we must compute the exponent in radix 2^32. We round it up
   * to the next multiple, and then shift left, in order to never lose
   * any information bits, adjusting the mantissa length if needed. */

  if (a_exp % 32 != 0)
    {
      int cnt = a_exp % 32;
      if (iexp < 0)
        cnt = 32 - cnt;

      ret = (iexp - cnt) / 32 + rl;
      if (uilsh (dst, dst, len, cnt))
        ++ret, ++rl;
    }
  else
    ret = iexp / 32 + rl;

  len = rl;
  return (ret);
}

double uitodbl (const limb_t *src, int len, int expo)
{
  int rexp = (expo - len) * LIMB_BITS;
  double ret = ldexp (*src, rexp);

  rexp += LIMB_BITS;
  for (int i = 1; i < len; ++i, rexp += LIMB_BITS)
    ret += ldexp (src[i], rexp);

  return (ret);
} 

KP_DECLS_END
